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SEMIDIRECT COMPUTATIONS FOR TRANSONIC FLOW 


Julie M. Swlsshelm and John J. Adamczyk 


National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 44135 


ABSTRACT 

A semldlrect method, driven by a Poisson solver, has been developed for 
In.vlscid transonic flw computations. It Is an extension of a recently Intro- 
duced algorithm for solving subsonic rotational flows. Shocks are captured by 
Implementing a form of artificial compressibility. Nonlsentroplc cases are 
computed using a shock tracking procedure coupled with the Ranklne-Hugonlot 
relationships. Results are presented for both subsonic and transonic flows. 
For the test geometry, an unstaggered cascade of 20 percent thick circular arc 
airfoils at zero angle of attack, shocks are crisply resolved In supercritical 
situations and the algorithm converges rapidly. In addition, the convergence 
rate appears to be nearly Independent of the entropy and vortlclty production 
at the shock. 


INTRODUCTION 

Semldlrect algorithms are constructed by Incorporating a direct solver In- 
to an Iterative procedure. For transonic flows, one of the first publications 
on this topic Is Lomax and Martin^). There, the governing equation was the 
small disturbance transonic flow equation. Jameson^ 2 ) developed a semldlrect 
solver for the potential flow equation. It was based on a direct solver for a 
Poisson equation. For subcrltlcal flows this algorithm exhibited good conver- 
gence; however, when the flow became supercritical convergence of the scheme 
could not be achieved without the Injection of line relaxation. 

More recently, a general methodology for constructing such schemes for the 
equations of fluid motion was advanced ty Mart1n( 3 ). The applications In his 
paper, however, were restricted to Incompressible and subsonic two-dimensional 
potential flows. Chang and Adamczyk^ 4 ) have succeeded In extending the scheme 
suggested by Martin for solving two-dimensional potential flows to three- 
dimensional Invlscld subsonic rotational flows. 

In the present paper, Chang and Adamczyk' s method Is extended to the com- 
putation of supercritical flows. The satisfaction of both the continuity equa- 
tion and the vortlclty-veloclty kinematic relation Is achieved by recursively 
applying a direct solver to two Poisson equations In computational space. For 
Isentroplc flows, shocks are captured by Implementing the artificial compressi- 
bility technique presented by Hafez, South, and Murman^ 5 ). The algorithm Is 
also capable of treating flows with nonlsentroplc shocks by adding a shock- 
tracking operator coupled with the Ranklne-Hugonlot relations. 

Numerical experiments conducted during the development of the present al- 
gorithm showed that the asymptotic convergence history was governed by a real 
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eigenvalue. This? Implied that the algorithm could be accelerated by a power 
method as suggested by the work of Hafez and Cheng(&). 

Convergence histories are presented for the unaccelerated and accelerated 
form of the algorithm. These results Include Invlscld subcrltlcal and super- 
critical potential and nonpotentlal flows. In addition the Mach number x ;r1- 
butlons associated with a cascade of unstaggered 20 percent thick cl ret ... arc 
airfoils operating at subcrltlcal an' supercritical conditions are presented. 


ALGORITHM FORMULATION 
Steady Potential Flow 

In steady potential flow calculations, flows are assumed to be Isentroplc 
and Isoenergetlc, Implying Irrotatlonallty. This assumption has proven useful 
In approximating many real flow situations. 

The differential equations governing potential flow are: 
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0) 

f x iJ = 0 

(2) 


where p Is the density and iJ Is the velocity vector. Equation (1) repre- 
sents the continuity equation and equation (2) the Irrotatlonallty condition. 

A thermodynamic relationship which expresses density In terms of velocity Is: 
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where P 0 , T 0 , R, C p , y and U 2 are total pressure, total temperature, the 
Ideal gas constant, the specific heat capacity at constant pressure, the ratio 
of specific heat capacities, and the speed squared, respectively. 


The solution Drocedure for subcrltlcal flows Is defined by the following 
set of equations^ 4 ) , written In tensor form: 
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where F k , U k , g kl , g^, g and x k are the contravarlant mass flux vector 
density, the covariant velocity component* the contravarlant metric tensor, the 
covariant metric tensor, the determinant of the covariant metric tensor, and 
the computational coordinates, respectively. The use of repeated Indices 
denotes the standard summation convention of tensor analysis. For three- 
dimensional flows the free Index k takes on the value of 1, 2, and 3, while 
for two-dimensional flows Its values are restricted to 1 and 2. The variables 
<p and a are Interpreted as scalar correction quantities. Equations (4). 

(5), and (6) satisfy the continuity equation (1), generating the vector F k 
which Is divergence free. Likewise, equations (7) and (8) generate a vector 
U k which satisfies the Irrotatlonallty condition (eq. 2). This algorithm 
can be thought of as a two-step recursive scheme In which the first step gener- 
ates a solution to the continuity equation while leaving the Irrotatlonallty 
condition unsatisfied. The, second step corrects the Intermediate velocity field 
so as to satisfy the Irrotatlonallty condition. Upon convergence, the Inter- 
mediate vector U k approaches uj^ . The Implementation of the above algo- 
rithm requires the solution of two Poisson equations (l.e., eqs. (4) and (7)) 

In computational space for each Iteration cycle. These equations are readily 
solved by direct solvers or fast Poisson solvers. 

The Iteration procedure Is Initialized by assuming a uniform covariant 
velocity field and constant density, which Is sufficient to determine the quan- 
tity F k ( n_1 ). Applying a direct solver to equation (4), the scalar <p Is 
computed for the entire flowfleld. Intermediate values of the mass flux vec- 

— 4 / rst 

tor, F , and the velocity vector, U^, are calculated by equations (5) and (6). 
The direct solver Is then Invoked a second time to find the a field which 
satisfies equation (7). The corrected velocity vector Is then determined 
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from equation (8). At this point the density Is updated using equation (9) and 

the current velocity. The mass flux vector F^ n ^ Is recomputed according to 
equation (10) and returned to equation (4), which Initiates the next cycle of 
the Iteration scheme. 

Global conservation of mass Is maintained by staggering the variables on 
the computational mesh. As shown In figure 1, the quantities <p and a are 
placed at the (1J) th mesh points, while all other variables (such as p, U|<, 
and s) are calculated at the center of each mesh cell side. For equations (5) 
and (8), In which F k and are computed using the correction variables 
and c, the differencing scheme Is : 
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where A k and B represent F k and q> for equation (5) and, similarly, 
and a for equation (8). 


In transonic flow computations, artificial compressibility Is employed as 
a shock-capturing procedure. Density Is modified by an upwlnded differencing 
procedure, according to the equation: 


where 



+ avl U 
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and where a Is a constant of order 1. The variable M. , Is the Mach number 
th 

at the (1,j) mesh point. The quantity p ^ j replaces the density In the 
calculation of the flux and the velocity vector. It should be noted that In 
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adapting this form of artificial compressibility we have assumed that the 
streamwlse grid lines (x 2 * constant) nearly coincide with the streamlines. A 
similar procedure Is applied to the speed squared, q 2 , to damp a Mach number 
overshoot that occurs Immediately upstream of the 

-2 2 . , /_2 

q i.j " q u V’i-i.j 

where 

2 kl M ,, 
q = g U k U v 


p 

The coefficients v and a' are the same as defined previously, q Is the 

—2 2 

speed squared, and q replaces q In equation (9) for the density calcula- 
tion. 


shock: 

- q ? J < 16 > 


The measures of convergence, or the "residuals" of the computation, are 

taken to be the magnitudes of the quantities v • (p3) and v • (U - 
Both of these are driven to zero by the Iteration procedure. Upon convergence 
the vector F* will be divergence free and Irrotatlonal as prescribed 
by the governing equations. 


Nonlsentroplc Flow 

The potential flow formulation Is extended to permit the computation of 
two-dimensional Isoenergetlc flows with nonlsentroplc shocks. This Is accom- 
plished by Introducing a shock-tracking operator coupled with the Ranklne- 
Hugonlot relationships and the entropy transport equation. 

The governing equations for Invlscld, Isoenergetlc flows with nonlsen- 
troplc shocks are: 
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and the thermodynamic relationship 
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The n onlsentroplc Iteration scheme Is basically the same as that for Isentroplc 
flow.,, with some modifications to be mentioned below. Density, as defined by 
equation (15), Is reevaluated for the nonlsentroplc case according to the 
expression: 


-( S-S )/C 
a — ' o' v 


( 21 ) 


where C v Is the specific heat capacity at constant volume, and (S-S 0 ), the 
entropy rise across the shock, Is obtained from the Ranklne-Hugonlot relation- 
ship: 
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which, when written solely In terms of M n , the Mach number normal to the 
shock, becomes 
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The switching operator defined earlier Is used to track the shock 

location and Its geometry. If we consider the flow to have a uniform entropy 
distribution across the Inlet, the shock provides the only mechanism by which 
entropy can be Introduced Into the field. So Introduced, the entropy (accord- 
ing to eq. (19)) will remain constant along streamlines downstream of the 
shock. Since we have assumed the grid line x 2 ■ constant nearly approxi- 
mates a streamline, the entropy will be held constant along each streamwlse 
grid line downstream of the shock. Once the entropy field Is established, the 
density Is reevaluated according to equation (21). The vortlclty field gener- 
ated by the shock may be computed from the following equation^ 7 ): 


to 

P 



(24) 


where <*>, p, and x 2 are the vortlclty, the pressure (which Is related to den- 
sity and entropy by the Ideal gas equation of state), and the transverse compu- 
tational coordinate, respectively. The vortlclty field Is used to construct a 

covariant vector whose only component Is (l.e., = 0) according to 


A< n) = J dx 2 (25) 
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One should note that the curl of the constructed covariant vector Is equal 

to u. By adding the quantity (a^ - to ulj n ^ at the end of each 

iteration cycle the resulting velocity field upon convergence will satisfy 
equation (18) . 

For subsonic flows the stability bounds of the current algorithm were de- 
fined by Chang and Adamc.zyk( 4 ) . Based on this analysis and a series of com- 
putational experiments It appeared that the asymptotic convergence rate of the 
current algorithm Is determined by a real eigenvalue. This Implies that the 
algorithm can be accelerated by a simple application of the power method as 
outlined by Hafez and Cheng( 6 ). To analyze this acceleration procedure let x 
denote the dominant eigenvalue, and assume that the asymptotic convergence his- 
tory of the current algorithm Is governed by the equation: 

i k+1 * Xc k (26) 


where Is the error vector at the end of the k*h Iteration cycle. For 
the present analysis Is defined as: 

c k = B k - B (27) 


where b‘ k represents the value of the unknown vector B at the end of the 
k th Iteration cycle. The value of X may be estimated from the equation: 

X = s|B k+1 - B k l / E ls k - B k_1 i (28) 


where e denotes the summation over all components of B. With X known the 
limit of B k Is estimated by means of the equation: 


o B k B k+1 - _B k 
B « B + 1 _ ^ 


(29) 


Numerical results will show that this simple acceleration procedure Is 
most effective In accelerating the current Iteration procedure. 


RESULTS 

The algorithm has been tested for both subsonic and transonic flows 
through a two-dimensional straight channel with a 10 percent thick circular arc 
airfoil mounted on one wall. This simulates a cascade of 20 percent thick un- 
staggered circular arc airfoils at zero angle of attack. The computational 
grid, shown In figure 2, with 60 mesh Intervals In the lengthwise direction and 
16 Intervals In the transverse direction, Is orthogonal at the boundaries, as 
required by the present treatment of boundary conditions In the Poisson solver. 
The Poisson solver used In this study was constructed using block trldlagonal 
Inversion. The boundary conditions required the specification of the mass flow 
at the Inlet and the flow angle at the exit. Along the walls of the channel 
flow tangency was required. 
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Potential Flow Results 


The results for the Isentroplc flow computations are shown In figures 3 
and 4. For the subcrltlcal case with an Inlet Mach number of 0.5, the Isomach 
distribution In the channel and the Mach number along the upper and lower walls 
are shown In figures 3(a) and (b), respectively. B»;th figures Illustrate the 
symmetry of the computed flowfleld about midchord of the airfoil. 

The convergence measures, defined earlier as the divergences of the vec- 
tors (p3) and 0 - at each grid point, are shown In figure 3(c), which 

shows the average value of v • (plJ) at each Iteration for the unaccelerated 
and accelerated versions of the algorithm. The convergence history of 

v • (u - exhibits similar behavior, meaning that the algorithm conver- 

ges to the correct limit. The residual of the accelerated computation has L >en 
reduced ten orders of magnitude In about 50 Iterations. 

For the supercritical Isentroplc flow case, an Inlet Mach number of 0.675 
causes the formation of a shock In the channel at about 75 percent of the air- 
foil chord as illustrated by the Isomach contour distribution of figure 4(a). 
The shock Is resolved between two consecutive grid points, as one can see from 
the plot of surface Mach number In figure 4(b). The unaccelerated and accel- 
erated convergence histories of the divergence of the mass flux vector are 
plotted In figure 4(c). The accelerated version Is reduced ten orders of mag- 
nitude In 90 iterations. 


Nonlsentroplc Flow Results 

For the nonlsentroplc computations, a uniform entropy distribution across 
the Inlet Is assumed. Hence subcrltlcal computational results are Identical to 
those of the isentroplc case displayed earlier. In supercritical flows the 
shock Is tracked by monitoring the switching operator for large variations. By 
assuming the shock orientation to be normal to the streamwlse gridlines, the 
resulting entropy rise across the shock Is calculated. The flowfleld result- 
ing from an Inlet Mach number of 0.675 Is Illustrated by the Isomach contours 
In figure 5(a). The corresponding distribution of Mach number along the walls 
Is plotted In figure 5(b). One can see from these figures that the Inclusion 
of entropy and vortlclty In the Iterative scheme alters the shock strength 
slightly from that of potential flow, which In turn affects the flow down- 
stream. The convergence histories of accelerated and unaccelerated computa- 
tions, shown In figure 5(c), show that the residual is reduced ten orders of 
magnitude In 110 Iterations. The nonlsentroplc flow convergence rate Is 
slightly slower than that for potential flow. 


CONCLUSIONS 

A semldlrect method for applications to transonic flows with shocks has 
been presented. 

The ability of the algorithm to resolve flows modelled by the potential 
equation and the Euler equations for flows with constant total enthalpy Is 
demonstrated for subcrltlcal as well as shocked supercritical flows. 
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The applicability of the scheme to these cases has been sustantlated by 
the ror»Muxatiunal results presented In this paper. Shocks are resolved sharply 
anr: the resulting flowflelds compare favorably to those obtained by others for 
tne same problem. 

For the standard, or unaccelerated, version of the algorithm, converged 
solutions were attained In 130 to 210 Iterations. By accelerating the solver 
using a procedure which annihilates the dominant eigenvalue, the number of 
cycles Is reduced to between 50 and 170 for the same convergence criteria. 
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(a) Isomachs. 



AXIAL DISTANCE 


(b) Surface Mach number distribution. 



(c) Convergence histories. 

Figure 3. - Subcritical flow case, inlet Mach 
number - 0.5. 












